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We apply the coupled cluster method (CCM) to the Hamiltonian version of the 



latticised 0(4) non-linear sigma model. The method, which was initially developed 

OV 

for the accurate description of quantum many-body systems, gives rise to two distinct 



approximation schemes. These approaches are compared with each other as well as 



with some other Hamiltonian approaches. Our study of both the ground state and 
collective excitations leads to indications of a possible chiral phase transition as the 
lattice spacing is varied. 
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1 Introduction 

The physics of particles interacting through the strong forces at low energies is dominated 
by pions (see, e.g., the discussion in Ref. H). In the hadronic sector the exchange of single 
and correlated pions between nucleons determines the long-range and intermediate-range 
behaviour of the nuclear force. The QCD Lagrangian, which is supposed to govern all 
strong interaction physics, has a chiral symmetry. Were this symmetry realised in nature, 
we would expect all hadrons to come in chiral pairs in which each hadron has a partner of 
the same mass but opposite parity. The fact that this does not actually occur is due to the 
effect of spontaneous breaking of the chiral symmetry which arises from the development 
of a condensate of quarks and gluons. The pions are interpreted as the Goldstone bosons 
which occur in such a scenario. If QCD were exactly chirally symmetric we would expect 
to find them at zero mass. The fact that their rest-mass is about 140 MeV/c 2 is due to the 
effect of the small masses carried by the quarks. This can be expressed in the Gell-Mann- 
Oakes-Renner relation 0, which relates the pion mass to the quark condensate and to the 



up and down quark masses, m u and rrid respectively, 

ml = - I — J {^ijj){m u + m d ). (1) 

Here (tpip) is the quark condensate in vacuum. 

In recent years a simple scaling argument by Brown and Rho ||] has sparked heated 
discussions about how mesons behave in a medium, where the quark condensate changes 
due to a partial restoration of chiral symmetry. The same question can be asked for finite 
temperature. One believes that for temperatures of around the order of the pion mass the 
chiral condensate will disappear and chiral symmetry will be restored Q. It is believed 
that in heavy ion collisions it is possible to inject enough energy to create such a phase. 
It is not yet clear whether such a plasma can reach thermal equilibrium, nor whether the 
time scale is long enough for the chiral phase to develop. 

All these issues have renewed interest in the dynamics of systems of pions. At the same 
time a revival has occurred in so-called baryon chiral perturbation theory ||], where the 
most general effective Lagrangian compatible with chiral symmetry is constructed, and all 
constants are fitted to experiments. This approach, which has met great success in the 
meson sector [Q], is now actively being pursued in the baryonic sector. 

As mentioned above, pions find their origin in the chiral symmetry, which can be related 
to the fundamental QCD Lagrangian. In the dynamical breaking of the symmetry the chiral 
partners of hadrons disappear and a set of three massless bosons, the Goldstone bosons, 
emerges. Although the mechanism is well understood, it is difficult to see explicitly how 
pions arise from the original QCD Lagrangian, especially since the pions are known to 
be composite particles. Therefore one often works with an effective Lagrangian which is 
compatible with the consequences of spontaneous symmetry breaking. The pion dynamics 
is somewhat restricted, since we have partial conservation of the axial current (PCAC), 
which limits the interactions we can build into our models. One model that satisfies these 
constraints is the 0(4) non-linear sigma model, which we shall study in this paper. 

The class of non-linear sigma models plays an illustrative role in many physical phe- 



nomena |], [j. In particular, the 0(4) version of the model can be used to describe the 
dynamics of pions. In that case the 0(4) symmetry can be broken explicitly, but the isospin 
0(3) sub-symmetry must be retained. Using a perturbative approach for the model in 2 
dimensions (1 space and 1 time) one finds that it exhibits both infrared and ultraviolet 
divergences 0. The latter are thought to be responsible for dimensional transmutation, 
the process whereby the dimensionless coupling constant acquires a dependence on a fun- 
damental length, such as the QCD scale in QCD. 

The short-distance behaviour has to be regularised, which is generally done by putting 
the model onto a lattice. For some finite lattice spacing one expects a phase transition 
between a system with essentially free rotators which interact weakly (at large lattice spac- 
ings), and a system of tightly bound rotators which exhibit only collective behaviour (at 
small lattice spacings) . In the former phase all excitations are massive and are characterised 
by the quantum numbers of the free rotator. In the latter case the excitation energy gap 
disappears. However, for two dimensions (one space and one time) the situation is different. 
From exact results it is known that the phase transition disappears due to strong inter- 
actions among the Goldstone bosons H. The excitations remain massive, although their 
masses decrease with decreasing lattice spacing. This result is confirmed in dimensional 
regularisation H, as is the asymptotic freedom of the two-dimensional non- linear sigma 
model. 

Most results in 0(4) theory have been derived in a Euclidean framework ||. In that 
case there is also a clear correspondence with high-temperature spin models |10|| . In a 
Euclidean framework space and time are treated on an equal footing. Not only is the time 
coordinate discretised but it is also taken to be imaginary. Apart from the ground-state 
energy it is very hard to recover properties of the system, and it is not clear whether the 
Euclidean and the Minkowskian (or Hamiltonian) methods yield the same result in a non- 
perturbative framework, especially since the lattice regularisation is implemented in such 
different ways in the two methods. 

In a Hamiltonian framework time is taken to be real and continuous (as it is in nature) . 



One can construct quantum-mechanical (Schrodinger) states and calculate basic expecta- 
tion values. However, for gauge theories the gauge has to be fixed when constructing a 
Hamiltonian, and therefore gauge invariance might be destroyed in subsequent approxi- 
mations. Space is discretised on a lattice, and the lattice spacing is the only dimensional 
parameter in the model. The phase transition from the disoriented (symmetric) phase 
of free rotators to the oriented (spontaneously broken symmetry) phase of tightly bound 
rotators can be seen as the dimensional transmutation, which sets the scale at which the 
system starts to exhibit the basic continuum behaviour. 

As the non-linear sigma model is an effective field theory for pions these features should 
put pion dynamics in a broader context. Indeed, we see that in low-energy pion dynamics 
the interaction vanishes as the momenta go to zero, which follows rather naturally from the 
fact that the pions are the Goldstone bosons of the chiral symmetry and the interaction 
is dictated by the spontaneous symmetry-breaking mechanism. However, since the pion 
is a composite particle the connection should break down for higher momenta (preferably 
somewhere near the QCD scale). This makes it difficult to connect excitations in the large 
lattice spacing systems we are investigating to actual physical particles. The system is 
more closely related to the dynamics of rotators described above. 

Apart from its intrinsic physical interest the non-linear sigma model is an ideal testing 
ground for non-perturbative methods. Although the Lagrangian is fairly simple and the 
degrees of freedom are restricted, it still has the richness of a full physical theory, whose 
properties should be correctly described by a proper non-perturbative method for the study 
of lattice field theories. The structure of field theories on the lattice is such that it begs to 
be treated by many-body techniques which first allow the ground state to be approximated 



in a well-determined and controlled way. Recently, Chin [11] has performed such a many- 
body calculation for the model being studied in this paper, using a variational technique 
embedded in a more general correlated basis function (CBF) approach. By contrast with 
the techniques used here, the CBF method, which uses a generalized Jastrow-correlated 
wave function cannot easily be formulated as a set of algebraic equations. Instead, it needs 



to be solved by Monte Carlo techniques on a finite portion of the infinite lattice, just as 
in most lattice-gauge approaches. The limit to an infinite lattice needs to be taken at the 
end of the calculation, using heuristic or theoretical scaling arguments. 

Our approach to the problem is the so-called coupled cluster method (CCM). The CCM 
is now widely acknowledged to provide one of the most widely applicable and most powerful 
of all microscopic formulations of quantum many-body theory. In the many applications 
which have so far been made in such diverse fields as nuclear physics, quantum field theory, 



condensed matter physics, quantum magnetism, and quantum chemistry [12], the CCM has 
been found to provide results which are among the most accurate available. Indeed, it has 
been shown to give results comparable to those from large-scale quantum Monte Carlo 
calculations in those cases where the latter can be performed. The interested reader is 



referred to Ref. 12] for a survey of the method and an overview of its applications. 

The CCM is formulated in terms of cluster correlation functions. These correlation 
functions, although typically only correlating the particles (or spins, or fields) on clusters of 
lattice sites with a finite spatial extent, do nevertheless extend the correlated wave function 
over the infinite lattice. In the CCM, unlike in CBF and most other lattice calculations, 
the lattice may be treated as infinite from the outset, as we will see below. 

In particular, the CCM has recently been applied to simple U(l) lattice gauge field 



theory fill, 14], and has been shown to provide an interesting scheme for studying the 



properties of such models. In the current work we introduce a variant of the CCM that 
works with functions rather than with the more usual operators. This approach will be dis- 
cussed in considerable detail in the present paper, although we shall succeed in constructing 
an operatorial approach as well. 

The remainder of this paper is organised as follows. In Sec. [2] we briefly discuss the 
0(4) non-linear sigma model and the Hamiltonian we construct from it (with some details 
relegated to Appendix [A]). Then, in Sec. |3|, we shall succinctly introduce the traditional 
form of the CCM and, in some more detail, the functional form. In Sec. || we lay the 
theoretical framework for an investigation of the 0(4) non- linear sigma model using the 
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traditional CCM approach. The same is done in Sec. || for the functional method. In Sec. ^ 
we collect the results from the various schemes, and contrast these both among themselves 
and with the few other results available. Finally we give a summary and outlook in Sec. pi. 

2 The 0(4) non-linear sigma model 

The 0(4) non- linear sigma model is defined by the Lagrangian density 

C = ~Tr(^£/^£/t), (2) 

where U is an SU(2) matrix-valued field of 2 x 2 unitary matrices, which is often conve- 
niently represented by a sum of unit and Pauli matrices, 

U = n^I + it • n, (3) 

and the unitarity constrains the four-vector n = (n, n±) to have unit length. 

We shall study the model on a spatial grid which is a cubic lattice in three dimensions 
and with continuous time, in a Hamiltonian formalism. We replace the spatial derivatives 
by finite differences, with a lattice spacing a. Since the time-derivative part of the La- 
grangian describes a freely rotating four-vector constrained to the unit sphere, it is not 
implausible that the kinetic energy in the Hamiltonian form at each point in space can be 
reduced to the angular part of the Laplacian, as is shown in Appendix [A]. Another way to 
represent this result is 

The lattice Hamiltonian contains a trivial factor a . the volume of the unit lattice cell 
in D spatial dimensions. We implicitly assume that all expectation values scale with this 
volume. Here I a can be taken as either the left or right SU(2) generators, 

[led, Uj] = iT a U-Ai or [I ai , Uj] = iU-AjTa, (5) 
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and #ij = hi ■ hj is the angle between the two unit vectors h\ and hj describing U\ and U\. 
We use here the convention related to the 0(4) symmetry. Factors of a half might appear 
in other conventions related to the isospin symmetry (and see, e.g., Chin fill]). We note 
that these different conventions will also lead to corresponding alternative definitions of 
the kinetic energy. We ourselves choose the convention which yields the kinetic energy of 
a unit mass on the four-dimensional unit sphere. Relabelling A = -?, we have the classical 
spin model 

ff=i^I? + AX>-cos0ij), (6) 

i (U) 

where the sum over (ij) runs over all nearest-neighbour pairs, and counts each pair (or 

lattice link) only once. 

The Hamiltonian can thus be interpreted as that of freely rotating unit vectors (or 
quantum rotors), with a nearest-neighbour coupling that tends to align them. It is this 
competition between free rotation and an alignment force that makes these models inter- 
esting to study. It also shows that they are identical to classical spin systems on a lattice, 
which can be obtained as either the large- J limit of a finite- J lattice spin model, or in the 
calculation of the partition function of a classical lattice spin problem. 

3 Elements of coupled cluster theory 

We start with an O(N) Hamiltonian of the form 

F = ^K i + A^(l-n i -n j ). (7) 

Here K\ is the kinetic energy on the ./V-dimensional unit sphere at lattice point i. The 
eigenstates of K\ are the (hyper)spherical harmonics \Ja)- v 

tfil-Za^ejlJa),. (8) 

Here we have divided the quantum numbers into those (J) that determine the eigenvalues 
of K\, and those (a) that label the different degenerate eigenstates. 



There are now at least two distinct ways to formulate a coupled-cluster approach for 
such a classical spin model. Firstly, there is an operatorial approach, which models itself 
closely on the standard formulation of the CCM, albeit with rather artificial operators. 
Secondly, there also exists a functional approach, which was first applied to lattice QED 
problems. Although the latter functional approach appears to be more natural for the 
present study, both methods are investigated and compared below. 

3.1 The standard form of the CCM 

Coupled cluster theory is one of the mainstays of modern microscopic quantum many-body 
theory, and reviews of the method are available (and see, e.g., ref. Ifl^] ). We shall therefore 
only give a very short summary of the salient points of the method; for full details see 
Ref. [|l2| and references cited therein. 

The key feature of coupled cluster theory is its construction of a systematically improv- 
able approximation to the exact correlated ground state of an interacting system, by the 
action of the exponential of an operator S on some suitable model state or bare vacuum, 
|0), which may (but does not have to) be chosen as the exact ground state in some limit. 
No attempt is made to normalise this state against itself, but one rather constructs a bra 
state in a bi-orthogonal system such that the inner product of the bra and the ket states 
remains unity. One starts from a set of (generalized) multi-particle or multi-mode creation 
operators a n , where the corresponding hermitian conjugate annihilation operators destroy 
the vacuum, a n \0) = 0. We now write the correlated ground state as 




tf)=exp > S n ai |0). (9) 



This state satisfies an intermediate normalisation condition, (Ol^) = 1. In many-body 
applications, the exponential Ansatz of Eq. (||) guarantees proper size-extensivity and con- 
formity with the Goldstone linked-cluster theorem even when approximations are made. 
By contrast, the corresponding linear parametrisation which characterises the correspond- 
ing configuration-interaction method (or generalized many-body shell model) will not, in 
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general, obey either of these important properties when approximations are made. One 
also parameterises the corresponding correlated bra state as 

(*| = (0| (l + ^5 n aJexpf-^S n atJ ; (10) 

where the coefficients S n and S n are to be determined independently by an extremum 
requirement, which can be given the variational form 

S(I[S,S]) = 0, (11) 

with 

I[S,§\ = (*\H\V), (12) 

where I[S, S] is the functional with respect to the infinite set of coefficients S = {S n \ n = 
1, 2, • • • } and S = {S n \ n = 1, 2, • • • }. In practice, however, we need to truncate the sum 
over creation and annihilation operators, and a consistent truncation scheme is found by 
restricting the two sets of coefficients in the same manner. For this choice we also find that 
the Hellmann-Feynman theorem is satisfied, i.e. 

(V\d x H\y)=d x E, (13) 

where E is the ground-state energy, E = (\I/| H |\I/}. Finally, in a formulation of the theory 
where the coefficients S n and S n are time-dependent, we find that S n and S n are canonically 
conjugate variables. There is a deep connection between these three properties, and we 



cannot give up one without giving up something else as well [12]. In particular we note 
that although in untruncated, and hence exact, parametrisations of the form of Eqs. @ 
and ([n]), the hermiticity relation \ty) = \ty) /(^\^) would hold, when Eqs. (g) and ( |l0| ) are 
truncated this exact hermiticity is broken in general. Although Arponen |15| has discussed 
means to restore the hermiticity within the CCM formalism, we do not pursue this point 
further here. 
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3.2 The functional form of the CCM 

The functional form is similar to the method sketched above, and as far as we are aware 



has first been used in refs. [13, 14]. The key ingredient is the parametrisation of the 



ground-state ket wave function in the exponential form 

<{r}|tf> = exp[S(M)], (14) 

where \{r}) is some suitably chosen complete set of wave functions. Such an approach 
appears to be naturally suited to lattice Hamiltonians, where a complete set of functions 
can often be chosen to be the product of complete sets at each lattice site. Truncations, 
which in the operator form are performed on the number and type of creation operators, are 
now performed on the functional dependence of the function £({?-}), typically by expansion 
in a complete set of appropriate functions. Once we construct the mode of action of the 
Hamiltonian on functions of the variables {r} - typically a differential, but at worst an 
integro-differential, operator - we can proceed as before and use the respective bra-state 
parametrisation, 

<¥|{r}> = [l + S({r})]exp[-S({r})], (15) 

and a variational functional I[S, S] defined by 

I[S, S] = j d{r}d{r'}[l + 5({r})] exp [S({r})] ({r}\ H \{r'}) exp [S({r'})] . (16) 

One typically excludes a constant term from S and S (corresponding to the exclusion of 
the identity operator from the sum over the multi-configurational creation operators in 
the various terms in Eqs. @ and (|i~0|)). However, we now no longer have intermediate 
normalisation. Thus, the constant part of exp (<5({r})) is not one, since an arbitrary power 
of a function orthogonal to the constant function can contain a constant part. In this case 
we need to enquire what happens with the Hellmann-Feynman theorem and the canonical 
variables in the time-dependent formalism? The key to the answer lies in orthogonality. 
While this is provided for by Wick's theorem when using the operator form, in the func- 
tional form this requires an expansion in orthogonal functions, relative to (typically) a 
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product of integrations. This will be shown to be straightforwardly implementable in the 
lowest set of approximations (SUB2-re, where only n two-body correlations are retained) 
for the model under consideration here. Since many-body orthogonal polynomials are not 
straightforward to define, more complicated correlations involving more than pairwise cor- 
relations do not seem to be easy to formulate in the functional form of the theory. It may 
still be that for certain classes of lattice models a functional form is easier to use, and even 
though we lose some of the elegance of the standard CCM, we might nevertheless still be 
able to proceed using this method. 

4 The operatorial form of the CCM for classical spin models 

Let us first investigate the operatorial CCM approach when applied to classical spin models 
of the form given in Eq. (0) ■ We use the fact that the uncorrelated ground-state wave func- 
tion is given by the constant function (the product over the grid of the lowest eigenstates 
|0)j of the kinetic energy operator K{). We now define a set of excitation operators for the 
CCM by 

cL = l J «)i(°li ; J >°> ( 17 ) 

which are chosen such that their hermitian conjugates annihilate the vacuum. The quantum 
number J labels the energy eigenstates of a single rotor, and a labels the different states 
for a given energy. We shall use the notation 

i{n}>=ni^>> i{o»=ni°»>' < 18 ) 

i i 

where the states \n\) are the standard coordinate representation at lattice site i. We shall 
mainly be concerned with the SUB2 family of approximations, in which only correlations 
between pairs of spins are incorporated, and in which the coupled cluster operator can 
hence be written as 

S 2 = Y, &U = E si [4 ® Cjr] (0) , (19) 
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where I is now the label for the energy eigenstates in the relative angle between the spins 
at lattice sites i and j. The sum over [ij] is over all pairs of lattice sites, i and j, counting 
each pair once, by contrast with the sum (ij) which runs over the nearest-neighbour pairs 
only, and which we introduced previously in Eq. (|6|). We have suppressed the a labels, 
because these are implicitly summed over in the scalar coupling. Since the only scalar that 
we can construct from the two vectors n; and n, is their inner product, we find 

(^i,nj|52ij|0i,0j) = -=S 2 {hi-hj) = -==£ 2 (cos0 u ); (20) 

V"Af V"A r 

5 2ij = / dnidnj|ni,nj) ^==6 , 2 (cos6 , i j) (0i,0j|, (21) 

J V"iV 

where £In is the surface of the iV-dimensional unit sphere. The potential can be expressed 
in the same way. Since the potential only depends on the relative angle between nearest- 
neighbour spins we require only the scalar coupling of spins 



V- 



ij = / dnidrij\hi, fij) A (1 — cos6>ij) (rji,nj|. (22) 

We now deal with the operators S2 . The best way to calculate the expectation value 
of the variational functional 

I[S,S] = ({0}\(l + S)e- § He"\{0}), (23) 

is to use the usual nested commutator approach to the CCM, 

e- S He S = H+[H,S] + -[[#, S], S] + . . . , (24) 

which truncates at second order for the current Hamiltonian since the Hamiltonian in this 
form is a second order differential operator which can act only at one or two factors in the 
ket state. We then evaluate the ensuing expectation values by introducing complete states 
of intermediate states at both sides of H. 
We use completeness in the form 

\{m})= fdn({n})\{h})({h}\{m}), (25) 
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where J dfi({h}) describes the product, over the whole lattice, of integrals over the N- 
dimensional unit sphere appropriate to the general O(N) case. Relation (|25|) assumes the 
overlap 

(h\\rh\) = 5{h\ — rh\). (26) 

Using the fact that the operator S^, where we now and henceforth suppress the pairwise 
subscript 2, projects to the right onto the ground state for sites i, j, we find that the only 
parts of the commutators in Eq. (^) that contribute are these where all S operators are 
to the right of the Hamiltonian. 

Using the fact, as follows from normalisation, that 



MOO = J±. (27) 

where Q.^ is the surface area of the TV-dimensional unit sphere, i.e., O4 = 2tt 2 for 0(4), we 
get 

I[S, S} = J d^({n}) j d^({h'}) <{0}| 1 + S \{h}) ({n}\ H\{n'}) <{n'}|(l + S + U 2 ) |0) . 

(28) 
If N g is the number of lattice (or grid) points, we may now write 

({h}\ S |{0}) = (n N )- N ^ 2 J^ S^.ini ■ nj), (29) 

[ij] 
and 

(Ml S 2 1{0}) = (n N T N » /2 Y! S xi-j(«i • *})S Xt ,_ y (iH> ■ ny) (30) 

where the prime on the sum means that only totally disconnected diagrams are allowed, 
due to the projective properties of our creation operators. The notation Xi-j denotes an 
integer irrep label; it indicates that we assume that S is invariant under symmetries of the 
lattice, and we have to identify components at symmetry-equivalent lattice sites. Finally 

<{o}| 1 + s \{h}) = (n N y N ^ 2 1 1 + J2 V->i • «j) ) • (3i) 

14 



x o 

K 
S S V 

Figure 1: The basic building blocks for the diagrammatic analysis of the functional I. We 
denote S by a thin solid line, S by a thick solid line and the potential by a dashed line. 
A cross denotes the action of the kinetic energy operator. Open circles indicate that the 
coordinates at these lattice points are not integrated over. 

In a SUB2-type approximation (i.e., when the correlation operators S and S are two-body 
operators), we can use those parts of the Hamiltonian where only the relative variables, #ij 
defined by cos #ij = n; • n_j, are present 

H -+ E ^k) 9 ^ (sin2 w + E A ^ - cos *d- ( 32 ) 

where the kinetic operator K acts twice on each relative angle #ij: one from each lattice 
site, K\ and Ky This explains the factor of two in the kinetic energy terms in Eq. ( |32[ ) 
with respect to those in the original Hamiltonian defined in Eq. (||). 

The part of the integration measure over the 4-sphere dependent on 9 is sin 2 8 d6. 
Integrating over a single S n is the projection onto the model state, which is by construction 
zero. Using this fact, we get the sum of all possible connected diagrams, generated from 
the basic ones sketched in Fig. [j]. The key idea is to construct all diagrams for I[S, S] on 
the lattice that can be made from drawing any of the lines (S, S, V) between the lattice 
points, subject to the condition that the diagram is closed, no two functions S can connect 
at the same point, and the potential connects only nearest neighbours. Diagrams that do 
not contain a potential must contain one kinetic energy operator K, and they can contain 
at most one line S. 

In the so-called LSUB2 approximation we retain only nearest-neighbour pairwise cor- 
relations. In Fig. ^ we place the diagrams on a lattice. Since the ket-state CCM equations 
are obtained through variation of I[S, S] with respect to 5ij, we should not integrate over 
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Figure 2: The diagrams for the LSUB2 approximation, together with the non-linear equa- 
tion for S and the equation for the energy. 

the unit vector at these grid points. This is denoted by the open circles in the CCM equa- 
tion. For fixed position of the S line, the square diagram in / can be produced in various 
ways. We choose to pick one orientation, and give it a weight a} 1 labelling the number of 
equivalent diagrams. The square diagram is the only non-local diagram in /, i.e., it involves 
nontrivial integrations over additional lattice sites, other than i and j from the variation 
with respect to S{j. 

Of course both the action and the energy are infinite for an infinite lattice. It is clear, in 
the LSUB2 approximation, that when the terms proportional to S vanish, which happens 
after minimisation, the energy can be expressed as a sum over all nearest-neighbour links 
of the three diagrams shown in the last line of Fig. g which do not contain a thick line. 
Using the standard assumption that each of these diagrams is independent of the (nearest- 
neighbour) sites i and j for the infinite lattice, due to translational invariance, we obtain 
an expression for the energy proportional to the number iVj of (nearest-neighbour) links. 
Therefore, in future, we shall always calculate the energy per nearest-neighbour link. 

A very natural basis for the expansion of all functions of 9 is the set of eigenstates of 
the kinetic energy (i.e., to that part of the Laplacian depending only on the azimuthal 
angle). As is well known, these are, in any number of dimensions, Gegenbauer polynomials 
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or hyper spherical functions. Since we shall use these polynomials in many places, we have 
dedicated Appendix || to a summary of their properties. For the 0(4) case studied here 
the appropriate Gegenbauer polynomials C*(cos#) are the Chebyshev polynomials of the 
second kind, U n (cos6). 

The general equation determining S, shown schematically in Fig. g in the case of the 
LSUB2 approximation, contains terms with zero, one, or two functions S (thin lines) 
linked by a potential. The integrations over the unit-vectors at the connecting points can 
be performed without any problem, including the weight factor l/f^ = l/(2ir 2 ) obtained 
from the normalisation of the ground state, as 

1 



Q dn{h 2 )S a {hi ■ h 2 ){h 2 ■ n 3 ) 

1 f j i- \o I n n^i(«2-"3) 



— / d[i(n 2 )S a (cos6i 2 ) 



2 
- cos(0i 3 )±S a ,i- (33) 

Here we have used the addition theorem for Gegenbauer polynomials, U n (x) = C^(x), as 
discussed in Appendix ||. We have defined the generalized Fourier coefficient as 

S al = - f sin 2 9d6 2 cos(6)S a {cos 9), (34) 

7T Jo 



where the factor two in front of cos 6 in the integrand of Eq. (34) originates from the use 
of the properly normalised hyperspherical polynomial, U\{x) = 2x. Since the integrand 
only depends on the azimuthal angle, the integration over the other angles simply yields a 
trivial constant factor, 

— / dfi(n 2 )f(ni-n 2 ) = — I siv? 6 12 d0 12 sin 4>d(j)dxf {cos 12 ) 

= - f sin 2 8 12 d9 12 f (cos 9 12 ) (35) 

as explained in more detail in Appendix [A]. 

By applying the addition theorem twice we find 

—2 / dfi(h 2 ) / dfi(h 3 ) S a (ni ■ h 2 )(h 2 ■ n 3 )S b (h 3 ■ n 4 ) = \(n 4 ■ hi)S a ,iS btl . (36) 
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We thus conclude that all two- and three-line diagrams containing one potential line are 
proportional to cos 9. 



4.1 The LSUB2 approximation 

We can now easily write down the non-linear equation for S± (i.e., the nearest-neighbour 
pairwise correlation function), 



Qi 



sin 



-do (sin 2 6d e Si) +A(l-cos0)(l + 5i) 



1 



\ai l cosOSf x , (37) 



where a^ 1 gives the number of times we can generate the box diagram for fixed position of 
the S line, and Q\ projects on all parts orthogonal to the constant function, 

(Qi/)(cos 9) = /(cos 6) - - I sin 2 9' d9'f( cos9'). (38) 

7T Jo 

We rewrite Eq. ( |3"7| ) in terms of Si = 1 + S\, to remove the inhomogeneity from the 
equation, as 



Q\ 



l 



sin 



-d e (sin 2 9d e Si) + A(l - cos 9)Si 



-Aqi cos 9 Si i . 
4 > 



(39) 



We now consistently drop the bar over S\, and have therefore to remember that henceforth 
we must satisfy the constraint 

2 
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sm z 9d9S 1 (cos9) = 1. 



(40) 



We can expand the solution to Eq. (£3y) in terms of Gegenbauer polynomials C*, which 
happen to coincide with the Chebyshev polynomials U n , 



Si(cos9) = y Si^ n U n (cos9). 



(41) 



These polynomials satisfy the equation 
1 



sin 



d e (sm 2 9d g U n {cos9)) = -n(n + 2)U n (cos9). 



(42) 
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Using the property 



xU n (x)U m (x) \J\ - x 2 dx = \8 n ,m±i, 



the equation for Si, namely Eq. (37), can be rewritten as the matrix equation 



/ 



-iA 





±A 3 + A 


-\\ 


-|A 


8 + A -\\ 



\ : 



where k is a convenient shorthand for 



\ 




f Si, 




( e \ 






£1,1 




K 






Sl,2 







) 




K'-- ) 




^ : / 



,11 



K = x^-st v 



(43) 



(44) 



(45) 



After its redefinition, we have the constraint that Si t o = 1. It appears that we have one 
more equation than we have unknowns, until we remember that the first equation (the 
one containing e) is not part of the CCM equations; rather, it has been added by hand for 
elegance. Of course we can always add an additional equation with one additional unknown 
(e). We then choose e such that 51,0 = 1> an d find that we just have enough unknowns to 
satisfy all equations. 

In order to solve these equations for the coefficient Si i and the energy we have to invert 
the matrix 

/ A -AA fl ... \ 



A 



|A 





+ A 


-|a o o 


iA 


8 + A -\\ 



V 



(46) 



/ 



Obviously we need only the the top corner coefficients of the inverse, all of which can be 
re-expressed in terms of a single one of the coefficients, 



(.l- 1 )^ = (A^) 21 = 2{A- 1 ] 



n 



2A~ 



(A- 1 ) 22 =4(A- 1 ) 11 -4A- 1 , 



(47) 
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which depend only on A. The top corner coefficient can be determined through various 
means, e.g., by the use of Cramer's rule: 

{A ~ 1)u = T > (48) 

where d n is the determinant of the matrix A with the first n rows and n columns removed, 
which satisfy the recursion relation: 

d n = [\ + n{n + 2)}d n+1 - -\ 2 d n+2 . (49) 

The solutions can be expressed elegantly in the form 

e = wo (A) — wi(A)k, (50) 

5i,i = Wi(A)(1 + 2k/A). (51) 



Here wq and uj\ are elementary functions of (A l ] 



n 






Substituting the explicit form of k from Eq. (|45|), we see from the equation for the energy 
in Fig. |2] that at self-consistency e = E/Ni, where iVj is the number of nearest-neighbour 
links on the lattice. We thus find 

£/iV; = wo(A) - gAwiCAJoPS?,!, (54) 

Si,! = Wi(A)(l + iai 1 ^). (55) 

Equation ( |55| ) can easily be solved, and shows a perfectly regular behaviour, 



"1 \ w i( A ) v wi ( A ) / 

The only acceptable solution is the one with the minus sign, since only in that case do we 
have zero energy at A = 0, where u>i(X) — ► 0. 
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I[S,S]= 




+2 




V^ +4 \ 



Figure 3: The diagrams for the SUB2-2 approximation in two space dimensions. Note that 
because of the redefinition of Si — ► 1 + Si, the two diagrams shown in Fig. g which contain 
the bare potential line, and both a potential line and an Si line, respectively, between 
nearest-neighbour lattice points, are incorporated in the remaining diagrams. We have 
also suppressed the lattice summation over the left-most link in each diagram, which was 
stated explicitly in Fig. g. 

4.2 The SUB2-n approximations 

We now discuss what happens if we add pairwise correlations beyond nearest neighbour. 
In a SUB2-n approximation we truncate at n pairwise correlations. Some of the relevant 
diagrams are sketched in Fig. [3| for two space dimensions in the SUB2-2 approximation 
where we include only two types of correlation functions, namely those between spins on 
nearest and next-nearest lattice points. The additional correlations {S n \ n > 1}, beyond 
nearest-neighbour correlations, do not contain direct products of the correlation with the 
potential. In S n , the label n refers to the irrep label introduced below Eq. (|29|). Therefore 
they will form closed loops with the potential in which only the lowest generalized Fourier 
coefficient survives, and all lead to differential equations of the form 

1 



sin 



-do (sin 2 6 dgS n ) = cos 8 n n ; n > 1, 



(57) 



with K n a quadratic polynomial in all of the S mj i's. Such an equation can easily be solved, 
and leads to 



Q n 

^n.l — -7T) 
D 



(58) 
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where all other Fourier-Gegenbauer coefficients are zero. Equation ( |58|) is thus really a 
quadratic equation which only involves the first Fourier coefficients. 

The only new diagrams contributing to the energy are the ones involving one S function 
and the kinetic energy operator (e.g., the last diagram in Fig. ||). Since all {S n \ n > 1} 
are proportional to the Gegenbauer polynomial of order 1, and are eigenfunctions of the 
kinetic energy, the derivative term is proportional to U\. When we integrate this over the 
relative angle we get a zero, due to the orthogonality of the polynomials. We thus conclude 
that the correlation functions {S n \ n > 1} do not directly contribute to the energy; rather 
they contribute only indirectly through the additional terms due to them in the equation 
for Si t i. 

In the end we obtain the coupled set of quadratic equations 

a 11 
e = wo (A) -a>i(A)A— — Sf )1; 

Si,i = wi(A) 1 + ^a^ + i^arV".,! > 

\ ' n n,m / 

Sn,l = ^«n%l + ^<C'Vu;''>l' (59) 

m m,l 

In order to understand the significance of the factors af m and a™ one should realise 

that we have assumed that our correlation functions are invariant under the symmetry 

transformations of the lattice (which would appear to be a sensible assumption for the 

ground state). We thus assign the same unique numerical label to all correlation functions 

of such a representation of the lattice symmetry group. The labelling function shall be 

denoted by Xn, and typically has the property that if |n| < |m|, Xn < Xm- The coefficients 

af 1 and af m then denote the statistical weights of the corresponding diagrams, i.e., the 

number of times we can fit such a diagram onto the lattice with fixed end-points. These can 

easily be determined from a combinatorial calculation. We assume that the functions are 

invariant under the symmetries of the cubic lattice, which leads to a great simplification. All 

vectors which are related to each other by space symmetry belong to the same irreducible 

representation, xi- We reserve Xi = 1 f° r the unit-length nearest-neighbour link. The 
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«? 


— ^7^, <W a _b"*Xb-c"lXc-a) 




a,b,c 


«r 


= 2_^ "«Xa-b"iXb-c"mXc-d"lXd 




a,b,c,d 



a-coefficients may then be expressed in terms of their multiplicities 

(60) 
(61) 

The sums over a, b, c, and d run over all lattice points, subject to the restriction a ^ c, 
a/d, and similarly for b. These restrictions originate in the projection part of the creation 
operators. 

4.3 Full SUB2 calculations 

For the case of quantum lattice spin models of relevance to (antiferro)magnetism it has 
been shown that the full SUB2 equations can often be solved by going to Fourier space 



16 1 . The corresponding equations to be solved in the present context are not very different 
from those studied in Ref. [jlq| . Unfortunately, however, the ui\{\) factor in the second line 
of Eq. (|59|) makes it impossible to obtain usable equations in reciprocal space, since the 



selfconsistency condition seems to be too complicated to allow a solution. 

4.4 Collective excitations 

An extremely important aspect of calculations of the type discussed here is the behaviour of 
excitations within the model. There are many different ways to perform such calculations, 
but one of the most appealing ones is through the study of small amplitude fluctuations 
around the ground state. This can be formulated in several equivalent ways, but as long as 
we are only interested in those excitations that have a wave function of similar structure 
to the vacuum it is easy to show that this corresponds to solving the linear random phase 
approximation (RPA) eigenvalue problem 

dy S 5,\L, M y) = E nUx). (62) 

8S{x)6S(y) 
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In the current problem this reduces to a simple matrix diagonalisation. The matrix takes 
on the block structure 



Xi Y X 

Y 2 X 2 



(63) 



with 



/ 3 + A -|A 



(Xi 



4A 8 + A -iA 



v ; 



+ 2Aa x S^iSj^Si^i. 



(64) 



' y 



In practice the matrix in this equation needs to be truncated, but for all values of A 
truncation at order 30 seems to be more than adequate. 
The matrix Y\ has one non-zero row: 



{Y 1 ) ln = \a n l /A + \arS m , 1 /A, 



(65) 



and Y<i has one non-zero column, 



(Y 2 ) nl = Xai/A + XaJTSm^/A. 



(66) 



Finally, the matrix X 2 has the simple structure 



{X: 



2)mn 



Aa^/4 + \a%S ltl /4 



(67) 



where we note that in Eqs. (pq)— (p7|) the summation over repeated indices is implied. Since 
the RPA matrix is not symmetric, it is not obvious that its eigenvalues are real, and in 
general they can be complex. 

5 The functional form of the CCM for classical spin models 

We have already seen that the operatorial method is driven by the potential term in 
the Hamiltonian, since most of the diagrams contain a potential line. Conversely, the 
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functional method is driven by the kinetic term. The differentiations in the kinetic term 
pull correlation functions down from the exponential functions in Eq. (jig ) and link them up 
at the lattice point at which the kinetic energy operator acts. The potential term simply 
commutes with the exponentials, and acts as an inhomogeneous term in the non-linear 
equations. 

Since the kinetic term contains two differentiations, it links only two correlation func- 
tions: 

K\f(h\-h m ,h\ -n k ) = 

- Am/(ni-n m ,ni-n k ) 

- Ak/(«r«m,nr%) 

- [n m -h k - (ni-h m )(h l -h k )]f( 1 > 1 \h l -h m ,h l -n ]i ) 1 (68) 

where D-^ is the azimuthal part of the kinetic energy acting on the relative angle #ij , 

Df {cosO) = [- sin 2 6d 2 cose + 3cos6d cose ]f{cos6), (69) 

and where f^' l \u,v) = d u d v f(u,v). 

5.1 The LSUB2 aproximation 

In the functional method we assume, in LSUB2 approximation, that 

<{n}|tt>=expf^S(0 u )J (70) 

We find that the CCM variational functional can now be written in the LSUB2 scheme as 

I(S,S)=Ni- f sin 2 6d6 Ul + S(0)) [-(S" + 2cot(6)S' + {S') 2 ) +V(6)] , (71) 

which only involves one-link quantities. We get no terms equivalent to the square diagram 
in the functional method, since the function S commutes with the potential. As before, 
all open diagrams disappear, even the one where the open line contains a derivative. For 
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example, one of those diagrams consists of two functions S' linked by an S, which form a 
closed triangle formed by linking two correlation functions, from the kinetic term acting 
at lattice point i, with Sn^x. It has a value proportional to 

J = I d hid h^d h\S'{ni ■ njJiS^nj ■ hi)[(hi ■ nkX^i ■ h\) — (h^ ■ hi)]S(hi ■ hu)- (72) 



We now choose a coordinate system where the axes are oriented such that 



n\ 








ra k 



COS $2 

sin #2 


V I 



ni 



COS #3 

sin 63 cos 03 

sin 6*3 sin ^3 

V J 



(73) 



The integration in Eq. ( |72|) then factorises into the product of two integrations over the 
relative angles, with one integration over the hi, 



J = — dfi(hi) / I I sin Oj sin (f)jdOjd(f)jdxj 

x S' {cos 6 2) S' {cos 8 3) {sin 6 2 sin 6 3 cos fa) S {cos 6 2) = 0, 



(74) 



which vanishes due to the trivial zero that arises from the integration over fa. 

For the special case of the LSUB2 approximation it is better to rewrite the expression 
explicitly in terms of the exponentiated quantities, 

2 



I(S,S) 



Ni- I sin 2 6d6 

IT 



(1 + S(0)}exp{-S(0)} 



--1--| sin 2 qL exp{S(0)} + V{6) exp{Sf((9)} 
sin a do au 



After redefinition of S and S, 

u = exp(S'), v = (1 + S) exp(— S) 
we have a standard variational functional 



I{u,v) = Ni- sin 2 dd0v {cos d) 

1 d . 2 „ d . ... .. E 



(75) 



(76) 



sin 



2 ede 



u(cos6). 



(77) 
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in which we have now introduced the energy E as a Lagrange multiplier since we have 
lost the intrinsic normalisation from the exponentials. The equations for u and v, which 
are obtained by variation with respect to v or u, respectively, are now identical (since the 
operator is self-adjoint) and we have to solve the equation 



-27^ sm ^ + A (l-cos0)-- 



^{(cos6O=O. (78) 



sin 2 0d# d6 v ' N L 

We substitute u(cos6>) = y(9 / 2) / sin(0), 9/2 = ip and obtain Mathieu's equation 



y" + (4 J -^ - A + l| - 2 [-2A] cos 2^ y = 0. 



(79) 



The only problem that remains is what boundary conditions we need to impose. Clearly 
the original function u depends on cos 9, and must be even in 9. Since we had to divide 
the function y by sin9 to get u, we find that y must be odd in 9. This immediately fixes 
the lowest eigenstate to occur at 

E = \-l + i b2 (-2\), (80) 



where 62 is the lowest odd characteristic value of Mathieu's equation with period ir [17]. 
Notice that we also get an estimate of the excitation energies for free, by replacing 62 by 
bin- This reflects the fact that in this case the result is fully variational, and could equally 
well have been obtained by performing a variational calculation, 

with the sum-over-links Ansatz for the wave function, 

<{n}|V)=J>^. (82) 



We note that the corresponding Jastrow-type variational calculation of Chin [11 1 differs 
from Eq. ( |82] ) only by the use of a product of exponentials instead of the sum of exponentials 
used here. This relation of the functional form of the CCM with a variational calculation 
is an accident of the LSUB2 approximation; it does not persist to other approximations in 
the SUB2-n or higher schemes. 
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5.2 The SUB2-n-m approximations 

As in the case of the operatorial form we also wish to incorporate additional longer-range 
correlations of two point nature. In this case we need to truncate both on the number n 
of such longer-range functions as above, but also on the number m of basis functions in 
which each such SUB2-n correlation function is expanded. The resulting approximation 
scheme is called the SUB2-n-m scheme. 

For the general SUB2-n-m case, where we assume functional dependence on an irre- 
ducible set of correlation functions, we cannot use the simplifying technique sketched above, 
but we shall have to deal with the full complexity of the problem. Since the reference state 
is an eigenstate of the sum of all the kinetic energy operators, thus corresponding to the 
weak-coupling limit, a natural set of functions is the set of Gegenbauer polynomials, as 
will be seen below. 

We expand S as 

S = JZ 5 xm(wi • n i+m ), (83) 

[im] 

where both sums run over all lattice sites. We will introduce the specific details and the 
exact specification of our basis functions later. For S we use the Gegenbauer polynomials 
we introduced earlier, 

oo 

S = ^ y^ s mjk U k (hi • n i+m ). (84) 

[im] fc=l 

If we insert these correlation functions, S and S, into our functional / we find: 



A T 

{0\U k (n o -n r 



c~ \-| — icv»u ■'■mj 



'm,A; 



sin 

\2\o' 



—rpe ( sin2 ed eS Xm ) 



-(l-(n -n m ) )S' Xm S' Xm 

-J2(n -h m - (n m / • n m )(n ■ «m') 5 x m _ m , 5 i m , 

m' 

+Si Xm vJ\0), (85) 
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The potential function Vj_j between nearest-neighbour sites can be expressed in Gegen- 
bauer polynomials as 

y ; _j = \U (hi -hi)- -AEMrai -hi). (86) 

Because the equations are translationally invariant we restrict ourselves to correlations 
between the origin, i = 0, and the other lattice sites. Since only the derivative of S appears 
in the expression it is more convenient to parametrise this derivative directly in terms 
of Gegenbauer polynomials, rather than by parametrising S itself. This may readily be 
archieved, since a truncation on S in terms of Gegenbauer polynomials implies a truncation 
of S' one order lower. Thus, we write 

oo 

S 'x™ = Y s Xm,kU k (n ■ n m ). (87) 

fe=0 

In this case the sum over the Gegenbauer polynomials starts at the constant function, Uq, 
since S' unlike S can contain a constant part. 

A tedious but straightforward calculation gives the following expressions for the varia- 
tional equations, 

° = /^^( I )[E5^x.J+i + (j' + 2)j Xn j-iM-(i) 

J 4 3=1 

1 oo oo 1 oo 3 

~ 2 Y Y s x»,l( s xn,3+l - s X r>,J+l+2) U j( x ) ~ 4 Y Y s x n ,3-i s x n ,i U j( x ) 

j=0 1=1 j=0 i=0 

i oo j-2 
+ 4 Y Y S Xn,i-2-iS Xn ,^ i (x) 
3=0 i=0 

1 \^Y^ ( U + 2 ) g Xn- m ,j-lgXm,j-l 

m j=0 v J\J ' J 

JjXn-mJ + lfXmJ+1 _ S Xn-m,j-l S Xm,j + l _ S Xn-m,j + l S Xm,j-l \ tt / \ 

(j + i)(i + 2) ^ (j + i)2 ( i + i)2 ; ^ ; 

+ \{U Q {x) -hj^x))], (88) 

where x = uq ■ n n , and where we have used the addition theorem for Gegenbauer poly- 
nomials. In order to solve these equations we also use the orthogonality relation for the 
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Gegenbauer polynomials: 

djj,{x 



n U i {x)U j {x) = 8 ij . (89) 

Note that the integration measure contains the normalisation constant from integration 
over the other (than the azimuthal) angles of the four-dimensional unit sphere. Of course 
in practical calculations the truncations of the infinite sums are important. These are 
performed by assuming that the coefficients s Xn j vanish for j larger than a cut-off value 
m. 

Once we have done that we can determine the energy per nearest-neighbour link from 
the expression, 

1 °° In 

E/N t = l/(2D)^[--^ Sxn)/ ( Sxn ,/- Sw+2 )- i ( Sxn ,o) 2 + A], (90) 

n 1=1 

where D represents the number of spatial dimensions. 

5.3 Collective excitations 

One way to calculate excitations is to use the so-called Feynman technique. Thus, the 
excitation energy for the state |\I/ e ) = X |\I/) can be evaluated by the double commutator 
of the Hamiltonian with the excitation operator, if we assume that X \fy) is an eigenstate 
of the Hamiltonian. In the CCM one often uses the same approach, even though it is 
no longer true that the bra and ket excited states are generated by the same excitation 
operator, X, and thus the approximations in this scheme are not easily controlled. For 
the case of the functional form of the CCM we use the operator that generates the first 
excited state in the limit of zero coupling constant. We can calculate the energy gap in 
this approximation analytically, knowing only the properties of the ground state. 

As usual in the functional form, the truncation is performed directly on the wave 
function rather than on the operators, and we assume the SUB2 form for the ground-state 
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correlation functions and the LSUB2 form for the excitation function, 
<{n}|* e > = Xexpl^S(cos0ij)J , 

<*e|{n}) = [l + ^ 1 5(cos9i j )jexp(-^ 1 S(cos0 ij )jX, (91) 

where the excitation function takes the LSUB2 form, 

X = Y j ni-h y (92) 

We now evaluate the expectation value of the second commutator of the excitation function 
X with the Hamiltonian, and find the usual result for the excitation energy, 

A£ = 4(/^(l+S)e-S[[ff,X],V) (/Mpild + S^y 1 ( 93 ) 

The numerator is the second derivative of the functional I[S, S] with respect to si : i. The 
denominator can be represented by a set of diagrams of one link terms, with cos 2 6^ and 
triangles containing two one-link terms, cos 9^ cos #kj j an d S. 

We can also calculate the excitations from a small-fluctuation RPA-type expansion 
around the ground state, similar to the method used in the operatorial approach. The best 
way to derive such an approach is by studying the small-amplitude limit of the classical 
equations of motions that can be derived from the quantum action functional 

A[S,S] = f dt f^fil(l + S)e- s [id t -H]e s 
Jo J Q 4 3 

= -i f 7 dt [ d ^fV § S - f dtI[S,S]. (94) 

Jo J Q A 9 Jo 

If we expand S and S in Gegenbauer polynomials, we obtain an action of standard canonical 
form, 

rT _ /-T 

dt / Sn,i$n,i / 



A[S,S] = -i / dtV]sn,iSn,i- dtl[s,s}. (95) 

^o _ , Jo 
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We thus need to change the expansion of S' in polynomials, to an expansion of S. This 
is a straightforward calculation using the recursion relation of the polynomials. In order 
to establish orthogonality we had to use the fact that S had SUB2 form, so that we could 
change to relative variables. For more complicated correlations this is no longer possible, 
and we lose the simple form of the time-dependent variational principle. 

By linearising the Hamiltonian equation of motion derived from A we see that we only 
need consider the symmetric second-order expansion of / around its stationary ground-state 
value, 



I[S,S] « I[S ,S } + 



5 2 I[S,S] 



5S5S 



. 5S5S 

So, So 



- /[S , Sol + »|_ )S ,M, (96) 

where we have introduced the matrix T (and see Eq. (|10S| ) in Appendix |B|) which imple- 
ments the change in basis when we expand S' or S in terms of Gegenbauer polynomials. 
The lowest eigenvalue of the functional varied with respect to S and S is the lowest exci- 
tation energy. 

6 Solution and Results 

6.1 Numerical methods 

Apart from standard matrix diagonalisations and inversions, the most challenging numer- 
ical problem encountered in the present work is tracing a solution to a parametric set of 
non-linear equations, where the coupling constant A is the running parameter. This prob- 
lem has been well studied by numerical analysts, since it occurs in areas of economics and 
engineering as well as in physics. There are various ways to solve such a problem, but one 
of the most robust and stable ones is based on a technique developed by Rheinboldt and 
collaborators Pq]. This technique combines a predictor-corrector method with a Newto- 
nian method to solve systems of non-linear equations, as one steps through parameter space 
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following the solution. We have used the implementation of this scheme in the PITCON 



program developed by Rheinboldt and Burkardt |L9]. 



6.2 LSUB2 Results 

For this approximation there is a large difference between the functional and operatorial 
methods. Even though the functional method has the same solution for any spatial di- 
mensionality, and is explicitly described by the solution to the Mathieu equation, this is 
not true for the operatorial method, where the result depends on the spatial dimension- 
ality through the statistical weight of the square diagram. Substituting a\ = 2 for two 
dimensions, and a\ = 4 for three dimensions, we find the results for the ground-state 
energy sketched in Figure [|. Both of these operatorial CCM solutions actually terminate 
and return. The termination points which occur for larger values than shown in Fig. |] 
are where us\ = 1/ ' y/a^ ■ ^ n one spatial dimension the non-local diagram does not exist, 
and hence there is no termination point in this case. However, in two and three spatial 
dimensions the termination point is at A = 115.2 and A = 16.3 respectively. By contrast, 
the functional form does not terminate and turn around. Rather, the solution exists for 
all A in this case, and hence no termination point exists. 

6.3 SUB2-n(-m) Results 

For the operatorial form we need to solve the non- linear Equations (^). Using the PITCON 
approach mentioned above, the equations can easily be solved once we determine the 
functions ujo(X) and uji(\), by truncating the set of equations where we keep only the 
lowest n coefficients Sk,i with k = 1, 2, • • • , n (and set the higher terms with k > n to zero). 
This approach is usually called the SUB2-n approximation. There is no real computational 
restriction on the value of n used in this approach, apart from the generation of the a- 
coefficients, which becomes very time-consuming for large values of the truncation index 
n. As we shall see below, convergence is attained before that becomes a real burden. 

In the functional method we need to truncate on the number of different S functions, 
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Figure 4: The energy per link in LSUB2 approximation, obtained in 2D (solid line) and 
3D (dashed line) from the CCM operatorial form, as well as by the solution to the Mathieu 
equation (dashed-dotted line) from the CCM functional form. 



corresponding to the set of distinct lattice vectors retained (i.e., those which differ under 
lattice symmetries), and for each distinct pairwise correlation function thereby retained 
we also have to truncate further on the number of basis functions in which they are ex- 
panded. Generally, the corresponding SUB2-n-m results depend more on the the number 
n of pairwise correlation components taken into account than on the number m of basis 
functions are used for the expansion of each component. Even close to the termination 
point all components behave smoothly and are not very large. They are well approximated 
with just three basis functions. 

We have studied the behaviour of the energy in 1, 2 and 3 space dimensions, and 
we present some relevant results in Figs. [5]-|7[ In each of these figures we compare the 
functional to the operatorial approach, both truncated at such high values of n and m 
that the results are converged on the scale shown in the diagrams. We find, following the 
path of the solution, that for all three cases the solutions are not single-valued in A. Each 
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Figure 5: The energy per link in the SUB2-ra (operatorial, solid line) and SUB2-n-m 
(functional, dashed line) approximations in 1 space dimension. Both schemes exhibit a 
termination point, which is denoted by the solid squares. The results are shown for suffi- 
ciently high values of n and m for convergence to have been attained on the scale of the 
figure. 



solution actually consist of two branches, that lie closer and closer together as we increase 
the order n. One of the two branches is "unphysical" in the sense that it does not connect 
directly with with the known solution for small coupling constants. This is also borne out 
by the excitation energies, as we see below. The mathematical reason for this behaviour is 
obvious. Thus, roots of non-linear equations cannot disappear; they have to collide with 
another root. Physically we interpret this behaviour as a sign of a break-down point of 
our approximation, and a possible indication of a phase transition at the corresponding 
termination point. The strength of this latter assumption is of course weakened by the fact 
that we find a phase transition here for the one-dimensional model where no such physical 
phase transitions occur |2(J. The fact that our putative phase transition occurs at smaller 
and smaller values of A as we increase the number of spatial dimensions may indicate that 
the occurrence of such a transition is more plausible in higher dimensions. 
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Figure 6: The energy per link in the SUB2-ra (operatorial, solid line) and SUB2-n-m 
(functional, dashed line) approximations in 2 space dimensions. Both schemes exhibit 
a termination point, which is denoted by the solid squares. The results are shown for 
sufficiently high values of n and m for convergence to have been attained on the scale of 
the figure. 



As we have noted earlier, the convergence to the limit point is fast in terms of the 
number m of basis functions required, probably indicating that the CCM is a good way to 
evaluate the ground-state properties. The number n of pairwise correlation components is 
the major factor in the convergence, as one would expect at a phase transition where the 
correlation length becomes infinite. In order to extrapolate the SUB2 critical value of A 
corresponding to the termination point, we perform a fit of the form 



Al = Aoo + c/L 2 , 



(97) 



where L is the distance between the correlated pairs which are furthest apart in the SUB2-ra 
calculation for a given value of the truncation index n. This asymptotic form is heuristic 
rather than having any theoretical basis. We start this extrapolation at L = 3, and make a 
standard linear fit to the numerical results, as indicated by the dashed lines in Figs. |-1C for 
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Figure 7: The energy per link in the SUB2-ra (operatorial, solid line) and SUB2-n-m 
(functional, dashed line) approximations in 3 space dimensions. Both schemes exhibit 
a termination point, which is denoted by the solid squares. The results are shown for 
sufficiently high values of n and m for convergence to have been attained on the scale of 
the figure. 

the operatorial cases. We contrast the results of the functional and operatorial approaches 
in Table [y. 

To our surprise, and somewhat to our dismay, we find that both the functional and 
operatorial methods predict a termination point, possibly indicative of a phase transition, 
in one dimension. This probably shows a defect in our approximation scheme, since no 
transition of this type is expected for such a one-dimensional system, and the exact results 
by Polyakov and Wiegmann || for the Euclidean theory also show no phase transition. The 
positions of the calculated termination points are surprisingly similar for the two approaches 
in all numbers of spatial dimensions, however. We also believe that both calculations have 
converged, so there is no obvious numerical explanation for the small differences. 

Of course we have not yet addressed the important issue of whether our results in- 
dicating a phase transition at finite A survive in higher-order approximations. It will be 
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Figure 8: The SUB2 termination point as a function of L 2 defined in Eq. ( |97| ) in 1 space 
dimension, using the CCM operatorial form. 

interesting to pursue this matter, especially for the ID case (where such a more sophisti- 
cated calculation to incorporate correlations between triplets or other n-body clusters with 
n > 2 may be more feasible than in higher spatial dimensionality). 



6.4 Excitations 

Excitation energies were evaluated by the linear response technique, for both the func- 
tional and operatorial CCM methods, as well as by Feynman's technique for the functional 



method only. The solutions for 1,2, and 3 dimensions are given in Figs. 11-13. In all three 
cases we find that the excitation energy obtained from the Feynman technique lies above 
that obtained from the corresponding linear response (RPA) technique. As we approach 
the termination points the excitation energies approach zero rapidly. However, since the 
numerical calculations become unstable at the the termination point, Feynman's technique 
only gives a very strong indication that at this point the system becomes gapless. 

The RPA linear response methods, where we expand around the ground-state solution, 



38 



c< 



H 1.3 




Figure 9: The SUB2 termination point as a function of L 2 defined in Eq. ( |97| ) in 2 space 
dimensions, using the CCM operatorial form. 



give slightly lower values for the excitation energy. In these calculations the excitation can 
be tracked all the way to the termination point, where the system becomes gapless. There 
is some doubt whether the extremely small difference between the point where the solution 
terminates and where the RPA frequency goes through zero is due to numerical artefacts 
or has some physical significance, indicating that only the full SUB2 calculation will be 
exactly gapless at the termination point. This might, a priori, be expected since only the 
full SUB2 calculation contains the arbitrarily long-range correlations needed precisely at a 
phase transition point. 

So what is the interpretation of the phase transitions that we have seen? We have 
investigated the eigenvectors for the mode that becomes (almost) gapless, but unfortunately 
the structure of these vectors, which tell us about the excitation operators to the new 
vacuum, do not seem to show any discernable regularity. This is probably related to the 
inability of the present method to describe the state beyond the putative phase transition. 

The study of the linear response also allows us to display most vividly the occurrence 
of two branches in the energy spectrum. In Fig. [14| we show the lowest excited state, 

39 




Figure 10: The SUB2 termination point as a function of L 2 defined in Eq. ( |97| ) in 3 space 
dimensions, using the CCM operatorial form. 

where we follow the solution of the ID SUB2-ra equations to the "termination point" and 
beyond. At this point the solution turns around, and we find a second solution for each 
value of A where the ground-state energy is extremely close to the one from the physical 
solution (so close that it is impossible to distinguish the differences in the figures). We 
can follow this solution until we return close to the point A = 0, where some of the CCM 
coefficients diverge. The excitation energy in this second branch is negative, a not very 
desirable situation. However, the fact that the second branch does not return to the known 
physical ground state at A = 0, and that the excitation energies are negative leads us to 
denote this branch as being unphysical, an artefact of the CCM. 



6.5 Comparison with other work 

There exist only two sets of results with which we can easily and directly compare our 



results, namely the calculation by Chin [11] using the CBF technique, and a calculation 



using the strong-coupling expansion by Sobelman [21]. 
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dimension 


Aoo 


functional 


operatorial 


1 
2 
3 


1.927 
0.827 
0.536 


2.088 
0.887 
0.552 



Table 1: The extrapolated CCM SUB2 termination points for the functional and operatorial 
methods. 



Sobelman has calculated the strong-coupling expansion in a Hamiltonian framework for 
a general O(N) system in an arbitrary number of spatial dimensions up to fourth order. 
If we use his results to calculate the coupling constant for which the mass gap disappears 
we find systematically larger values for the (2+1)- and the (3+l)-dimensional models, by 
a factor of about 3, than obtained by our results. In the third-order expansion there is a 
vanishing mass gap for (1+1) dimensions at the critical value A = 8, but this criticality 
disappears in the fourth order. There remains a dip in the excitation energy around A = 8 
where it crossed the zero axis before. Prom the exact results this seems to be the correct 
behaviour. However, the large change between third and fourth orders suggests that the 
expansion has not converged for the (1+1) model in this regime. For (2+1) and (3+1) 
dimensions the positions of the vanishing mass gap change somewhat, as one goes from 
third order to fourth order (and see Table ^). The convergence properties of the series are, 
however, completely unknown. 

Because of a different definition of the kinetic energy operator, the value obtained by 



Chin [11] for the critical coupling constant needs to be multiplied by a factor of four to 
agree with our definition of the model Hamiltonian given by Eq. (0). His results are for 
three dimensions only and he finds a critical coupling constant of A ~ 0.96 (after scaling) 
compared to our results Aoo = 0.54 from the SUB2 functional form and Aoo = 0.55 from 
the SUB2 operatorial operatorial. However, more strikingly, he finds a first-order phase 
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Figure 11: The excitation energy as a function of A in 1 space dimension in the CCM 
SUB2 approximations. The solid line is the RPA result from the operatorial approach; the 
dashed line the values obtained using the Feynman technique in the functional approach, 
and the dashed-dotted line the corresponding RPA result using the functional approach. 

transition, while all our results indicate a second-order phase transition. 

7 Conclusions and outlook 

Our results are an interesting albeit still limited approach to a nontrivial field theory 
which exhibits a phase transition. They show that we can construct two different CCM- 
like approaches for field theories, which are amenable to different sets of approximations. 
For similar truncations the two methods produce very similar results. 

In one space dimension we would not have expected a termination point or a phase 
transition from known exact results. The fact that in our approximations the termination 
point occurs at a large value of A shows that it is probably not stable against the correlations 
missing from our SUB2 calculations. The corresponding stability of our results in two and 
three spatial dimensions is certainly also of interest. The current calculations cannot 
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Figure 12: The excitation energy as a function of A in 2 space dimensions in the CCM 
SUB2 approximations. For the meaning of the lines see Fig. 11. 



address this issue, since they would have to be supplemented by complicated higher-order 
calculations in order to make stringent statements. This further raises the question whether 
one should not aim to apply the extended coupled cluster method (ECCM) |^2| to the 
present problem. The ECCM generally has been shown to be more powerful than the 
normal (NCCM) version of the CCM used here in describing states on both sides of a 
phase transition, and might give us an idea of its nature. We believe that this might be an 
interesting problem for a future study. 

Finally it would be appealing to look at real gauge field theories such as compact QED, 
or even QCD, in the light of our results. We plan to see whether the lessons we have 
learned in this work have any relevance for the methods of Refs. ]l3|, [l4|. We believe that 
armed with our arsenal of techniques we should be ready to tackle such problems, and even 
be able to proceed towards the inclusion of dynamical quarks in such calculations 
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Figure 13: The excitation energy as a function of A in 3 space dimensions in the CCM 
SUB2-n-m approximations at the highest practical truncation indices n and to. Note that 
(for three dimensions) convergence in the index n near the termination point is not fully 



attained. For the meaning of the lines see Fig. 11 
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A Quantisation 

The systematic quantisation of the non-linear sigma model is highly nontrivial. Using 
the components of the unit vector h as degrees of freedom, the constraint among them 



n 



1) leads to a second constraint from the time independence of the first constraint. 



The Dirac brackets, which restrict the quantised operators to the physical subspace, are 



not particularly difficult to deal with [23], but they lead to a form that does not readily 
lend itself to the kind of techniques we would like to apply. One can also avoid using 
constraints by using an exponential parametrisation, U = e tT ' e . However, this leads to 
even more complicated expressions. 
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Figure 14: The excitation energy as a function of A in 1 space dimension using the opera- 
torial form of the CCM in SUB2 approximation, showing both the physical and unphysical 
branches. 

We choose to use the unit vectors h as degrees of freedom, but avoid the use of con- 
straints by using an Euler angle parametrisation of these vectors, 



h = (sin 9 sin <j) sin %, sin 6 sin (j) cos \i srn & cos </>) cos $) > 



(98) 



with < 6 < it, < (p < tt, and < x < 27r. Since the Laplacian is separable in spherical 
coordinates, we find that the kinetic energy K at a single lattice point is just the half of the 
angular part of the Laplacian, which is proportional to the Casimir invariant of the relevant 
O(N) group. The kinetic energy can be derived from the generalized angular momentum 
operators 



Lij — —i(nid nj — Tijdrn), 



(99) 



where i,j = 1,2, ••• ,N. Specifically, it may can be expressed as (the O(N) Casimir 
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dimension 


Aoo 


third order 


fourth order 


1 
2 
3 


8. 
2.624 
1.637 


2.464 
1.492 



Table 2: The coupling constants for which the excitation energies vanishes, determined 
from Sobelman's strong coupling expansion. 

invariant or) the sum-of-squares of the angular momentum operators, 

N 



2K = Y^ Lfj = -de ~ (N - 2) cot 

i<j 



1 2 (iV-3)cot0 a 

n : <Jj, : — o - ^ <Jtb 



sin 



sin 



I a 2 + 

sin #sin <p 



(100) 



In many cases we shall only need the ^-dependent part, which can be rewritten in terms 
of cos 8 as 



(1 - cos' 9)d z cosfj + {N - 1) cos 9d cos6 . 



(101) 



B Gegenbauer polynomials 

The Gegenbauer polynomials, C"(x), are orthogonal polynomials on the domain [—1,1] 
with respect to the weight (1 — x 2 )" -1 ' 2 . The integration measure for the azimuthal angle 
on an TV-dimensional sphere is sin _ 9, which means that orthogonal functions on the N- 
sphere are C"(cos9) with a = (N — 2)/2. The Gegenbauer polynomials are eigenfunctions 
of the kinetic energy operator 

-[(1 - x 2 )d 2 x - (TV - l)x^]C^- 2 )/ 2 (x) =n(n + N- 2)C^ v ~ 2 )/ 2 (x). (102) 
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For TV < 4 there is a second, non-singular set of eigenfunctions of the relative angle 

-[(1 - x 2 )d 2 x -(N- l)x0 x ](l - x 2 )^- N ^ 2 Ci 4 - N ^ 2 (x) = (n + l)(n - JV + 3)C< 4 -^/ 2 (x), 

(103) 

which coincides with the first set of solutions for N = 3, and leads to a set of sine functions 
for N = 2. For a general discussion of the properties of the Gegenbauer polynomials we 
refer to [24]. Here we shall just summarise a few key points. 



One relation that we shall use frequently in the paper is related to the addition theorem 
for Gegenbauer polynomials. This occurs in the integration over intermediate points in a 
chain of functions linking sets of lattice points. Only products of Gegenbauer polynomials 
of the same order contribute to the integral. This holds for arbitrary dimensions, and 
allows us to extend our methods to any O(N) theory. Using the addition theorem for 
Gegenbauer polynomials 

C° (cos 9\ cos 82 + sin 9\ sin $2 cos <fi) = 

r( 2 q - 1) " 2 2m T(n - m + l)(r(a + m)) 2 ^ 

r 2 (a) ^ r(n + 2a + m) [ ° m ' 

v ' m=0 ' 



x (sin 0!) 771 (sin 92) m C^ (cos 6!)C^ (cos 6 2 )C^ 2 (cos0), (104) 

and orthogonality of the Gegenbauer polynomials it follows that 

7 ~n^~ C * (ni ' nj)Q ^ ■ Hp) - T 2 (a)(a + k)T(a + 1/2) 5fcA (nj ' "^ (1 ° 5) 
When integrating over closed loops this general property implies that only "local terms" 
contribute, i.e., products of Gegenbauer polynomials of identical order. For N = 4 where 
a = 1 this reduces to 

— — U k (rii- nj)Ui(ni- n p ) — , (lUb) 

where C^(cos9) = U n (cos0) is the Chebyshev polynomial of the second kind. 
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The fact that the derivative of a polynomial is again a polynomial allows us to express 
the derivative of a Gegenbauer polynomial as a weighted sum of Gegenbauer polynomials 

[n/2] 

d x C%(x) = J2 2(n - 1 + a - 2i)C^ 2i _ 1 {x) . (107) 

i=0 

We shall apply this transformation for the case where n^O (we exclude a constant part 
from our correlation functions). In that case we can define a matrix T by 

[n/2] 

T nm = ^2 2 ( m + aO<W-2i-i, (108) 

8=0 



which in the subspace of interest is invertible. This linear operator T is used in Eq. (p&j 
in the RPA for the collective excited states. 
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